/***************************************************************************\
 * Copyright (C) by Keio University
 * BilateralFilter.cuh created in 01 2013.
 * Mail : fdesorbi@hvrl.ics.keio.ac.jp
 *
 * BilateralFilterTemplate.cuh is part of the HVRL Engine Library.
 *
 * The HVRL Engine Library is free software; you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or
 * (at your option) any later version.
 *
 * The HVRL Engine Library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 \***************************************************************************/

#ifndef HVRL_BILATERAL_FILTER_TEMPLATE_CUH__
#define HVRL_BILATERAL_FILTER_TEMPLATE_CUH__

#include <math.h>

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

#include <cuda.h>
#include <cuda_runtime.h>
#include <hvrl/gpgpu/cuda/Prototypes.hpp>

namespace hvrl {

namespace cuda {

/**
 * \brief This templated class apply a bilateral filter on an input image using CUDA
 *
 * Input images have to be grayscale images
 * The method initialize has to be called in order to setup the spatial part of the filter
 * \author Francois de Sorbier
 */
template<typename T>
class BilateralFiltering_CU {

public:
	/**
	 * \brief Constructor
	 * \author Francois de Sorbier
	 */
	BilateralFiltering_CU(void);

	/**
	 * \brief Destructor
	 * \author Francois de Sorbier
	 */
	~BilateralFiltering_CU(void);

	/**
	 * \brief Precomputes the spatial part of the bilateral filter kernel and few values of the range part.
	 *
	 * Description of the bilateral filter :
	 *
	 * \f$I_p = \frac{1}{W_p} \sum_{q \in S}{G_{\sigma_s}(\|p-q\|)G_{\sigma_r}(|I_p-I_q|)I_q}\f$
	 *
	 * with \f$W_p = \sum_{q \in S}{G_{\sigma_s}(\|p-q\|)G_{\sigma_r}(|I_p-I_q|)}\f$
	 *
	 * \f$G_{\sigma}(x) = \frac{1}{\sigma \sqrt{2\pi}}\exp({{-}{x^2}\over{2\sigma^2}})\f$
	 *
	 * \f$G_{\sigma_s}\f$ is named the space weight (distance from center) and \f$G_{\sigma_r}\f$ is the range weight (difference of intensity)
	 *
	 * \param sigma_s the standard deviation parameter of the Gaussian in the spatial part
	 * \param sigma_r the standard deviation parameter of the Gaussian in the range part
	 * \param radius the radius of the kernel
	 * \author Francois de Sorbier
	 */
	void initialize(const float& sigma_s, const float& sigma_r,
			const int& radius);

	/**
	 * \brief Apply the kernel filter onto the input image
	 *
	 *	The input image should be a grayscale image
	 *
	 * \param d_input the input grayscale image already on the device
	 * \param d_output the result of the filter on the device
	 * \param d_width the width of the image
	 * \param d_height the height of the image
	 * \param block_x the number of threads we want on the x axis of the cuda 2D block
	 * \param block_y the number of threads we want on the y axis of the cuda 2D block
	 * \author Francois de Sorbier
	 */
	bool apply(T *d_input, T *d_output, const unsigned int& d_width,
			const unsigned int& d_height, const unsigned int& block_x,
			const unsigned int& block_y);

private:

	float *d_spatial_kernel; /**< Constant size array for holding the spatial kernel located in the constant memory of the host */
	int d_radius; /**< Radius of the kernel located in the constant memory of the host */
	float d_coefficient_range; /**< Precomputed coefficient for the range filter located in the constant memory of the host */
	float d_divider_range; /**< Precomputed coefficient for the range filter divider located in the constant memory of the host */

};

template<typename T>
void bilateralFilterKernelGray(T *d_input, T *d_output, float* kernel,
		const int& radius, const float& range_coeff, const float& range_divider,
		const unsigned int& d_width, const unsigned int& d_height,
		const unsigned int& block_x, const unsigned int& block_y);

extern template void bilateralFilterKernelGray<unsigned char>(
		unsigned char *d_input, unsigned char *d_output, float* kernel,
		const int& radius, const float& range_coeff, const float& range_divider,
		const unsigned int& d_width, const unsigned int& d_height,
		const unsigned int& block_x, const unsigned int& block_y);
extern template void bilateralFilterKernelGray<unsigned short>(
		unsigned short *d_input, unsigned short *d_output, float* kernel,
		const int& radius, const float& range_coeff, const float& range_divider,
		const unsigned int& d_width, const unsigned int& d_height,
		const unsigned int& block_x, const unsigned int& block_y);
extern template void bilateralFilterKernelGray<float>(float *d_input,
		float *d_output, float* kernel, const int& radius,
		const float& range_coeff, const float& range_divider,
		const unsigned int& d_width, const unsigned int& d_height,
		const unsigned int& block_x, const unsigned int& block_y);
extern template void bilateralFilterKernelGray<int>(int *d_input, int *d_output,
		float* kernel, const int& radius, const float& range_coeff,
		const float& range_divider, const unsigned int& d_width,
		const unsigned int& d_height, const unsigned int& block_x,
		const unsigned int& block_y);

template<typename T>
BilateralFiltering_CU<T>::BilateralFiltering_CU(void) {
	this->d_spatial_kernel = 0;
	this->d_radius = 0;
	this->d_coefficient_range = 0.0f;
	this->d_divider_range = 0.0f;
}

template<typename T>
BilateralFiltering_CU<T>::~BilateralFiltering_CU(void) {
	this->d_radius = 0;
	this->d_coefficient_range = 0.0f;
	this->d_divider_range = 0.0f;
	if (this->d_spatial_kernel != 0) {
		checkStatus(cudaFree(this->d_spatial_kernel));
	}
}

template<typename T>
void BilateralFiltering_CU<T>::initialize(const float& sigma_s,
		const float& sigma_r, const int& radius) {

	float coefficient_space = 1.0f / (sqrtf(2.0f * M_PI) * sigma_s);
	float divider_space = -0.5f / powf(sigma_s, 2);
	this->d_coefficient_range = 1.0f / (sqrtf(2.0f * M_PI) * sigma_r);
	this->d_divider_range = -0.5f / powf(sigma_r, 2);
	this->d_radius = std::abs(radius);

	int width = this->d_radius + 1;
	float *kernel = new float[width * width];

	// Precompute the 2D spatial kernel
	for (int i = 0; i < width; ++i) {
		for (int j = 0; j < width; ++j) {
			kernel[i * width + j] = coefficient_space
					* exp(
							(pow(float(i), 2) + pow(float(j), 2))
									* divider_space);
		}
	}

	//Check if the data have to be allocated on the device, and allocate the memory on the device
	if (this->d_spatial_kernel != 0) {
		checkStatus(cudaFree(this->d_spatial_kernel));
	}
	checkStatus(
			cudaMalloc((void**) &(this->d_spatial_kernel),
					width * width * sizeof(float)));
	checkStatus(
			cudaMemcpy(this->d_spatial_kernel, kernel,
					width * width * sizeof(float), cudaMemcpyHostToDevice));

	delete[] kernel;

}

template<typename T>
bool BilateralFiltering_CU<T>::apply(T *d_input, T *d_output,
		const unsigned int& d_width, const unsigned int& d_height,
		const unsigned int& block_x, const unsigned int& block_y) {

	if (this->d_spatial_kernel == 0)
		return false;
	bilateralFilterKernelGray(d_input, d_output, this->d_spatial_kernel,
			this->d_radius, this->d_coefficient_range, this->d_divider_range,
			d_width, d_height, block_x, block_y);
	return true;
}

} // End of the namespace gpgpu

} //end of the namespace hvrl

#endif
